<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>A step-by-step guide to parsimoneous enzyme usage Flux Balance Analysis - pFBA</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S2 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S3 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S4 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }
.S5 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S6 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>A step-by-step guide to parsimoneous enzyme usage Flux Balance Analysis - pFBA</span></h1><div  class = 'S1'><span style=' font-weight: bold;'>Author(s): Francisco José Pardo Palacios, Ines Thiele, LCSB, University of Luxembourg.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Reviewer(s): Sebastián Mendoza, Center for Mathematical Modeling, University of Chile. Catherine Clancy, LCSB, University of Luxembourg.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Lin Wang (Costas D. Maranas Lab), Joshua Chan (Costas D. Maranas Lab), Chiam Yu Ng (Costas D. Maranas Lab)</span></div><h2  class = 'S2'><span>INTRODUCTION</span></h2><div  class = 'S1'><span>This tutorial shows how the parsimoneous enzyme usage Flux Balance Analysis (pFBA), as described in Lewis et al.</span><span texencoding="$^1" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>, has been implemented in The COBRA Toolbox as the function </span><span style=' font-family: monospace;'>pFBA()</span><span>.</span></div><div  class = 'S1'><span>The main aim of the tutorial is to explain how the calculations are carried out in order to understand the pFBA analysis, and to be able to classify, under certain conditions, the genes of a model as: essential, pFBA optima, Enzymatically Less Efficient (ELE), Metabolically Less Efficient (MLE) or pFBA no-flux genes (Figure 1). </span></div><ol  class = 'S3'><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>Essential genes:</span><span> metabolic genes necessary for growth in the given media.</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>pFBA optima: </span><span>non-essential genes contributing to the optimal growth rate and minimum gene-associated flux.</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>Enzymatically less efficient (ELE): </span><span>genes requiring more flux through enzymatic steps than alternative pathways that meet the same predicted growth rate.</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>Metabolically less efficient (MLE):</span><span> genes requiring a growth rate reduction if used.</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>pFBA no-flux:</span><span> genes that are unable to carry flux in the experimental conditions.</span></li></ol><div  class = 'S5'><img class = "imageNode" src = "" width = "479" height = "481" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>                                                                                 Figure 1: Gene/enzyme classification scheme used by pFBA</span></div><div  class = 'S1'><span>This tutorial will use the</span><span style=' font-style: italic;'> E. coli core </span><span>reconstruction</span><span texencoding="$^2" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span> as the model of choice, and will be called herein as </span><span style=' font-family: monospace;'>modelEcore</span><span>. The results obtained could then be compared to data from evolved E. coli and observe if the </span><span style=' font-family: monospace;'>modelEcore</span><span> can predict its evolution. In order to investigate this, all the steps described in the pFBA flowchart (Figure 2) should be followed, and are demonstrated in this tutorial. </span></div><div  class = 'S5'><img class = "imageNode" src = "" width = "506" height = "587" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>                                                                                                          Figure 2: pFBA flowchart</span></div><h2  class = 'S2'><span>EQUIPMENT SETUP</span></h2><h2  class = 'S6'><span style=' font-weight: bold;'>Initialize the COBRA Toolbox.</span></h2><div  class = 'S1'><span>If necessary, initialize The Cobra Toolbox using the </span><span style=' font-family: monospace;'>initCobraToolbox</span><span> function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div></div></div><h2  class = 'S6'><span style=' font-weight: bold;'>Setting the </span><span>optimization</span><span style=' font-weight: bold;'> solver.</span></h2><div  class = 'S1'><span>This tutoria</span><span>l will be ru</span><span>n </span><span>with a </span><span style=' font-family: monospace;'>'glpk</span><span style=' font-family: monospace;'>'</span><span> package, which is a linear programming ('</span><span style=' font-family: monospace;'>LP'</span><span>) solver. The </span><span style=' font-family: monospace;'>'glpk</span><span style=' font-family: monospace;'>'</span><span> package does not require additional instalation and configuration.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >solverName = </span><span style="color: rgb(170, 4, 249);">'glpk'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >solverType = </span><span style="color: rgb(170, 4, 249);">'LP'</span><span >; </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >changeCobraSolver(solverName, solverType);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% changeCobraSolver('ibm_cplex');</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% changeCobraSolver('gurobi');</span></span></div></div></div><div  class = 'S11'><span>However, for the analysis of larger models, such as Recon 2.04</span><span texencoding="$^3" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>, it is not recommended to use the </span><span style=' font-family: monospace;'>'glpk</span><span style=' font-family: monospace;'>'</span><span> package but rather an industrial strength solver, such as the </span><span style=' font-family: monospace;'>'gurobi' or 'ibm_cplex'</span><span> package.</span></div><div  class = 'S1'><span>A solver package may offer different types of optimization programmes to solve a problem. The above example used a LP optimization, other types of optimization programmes include; mixed-integer linear programming ('</span><span style=' font-family: monospace;'>MILP</span><span>'), quadratic programming ('</span><span style=' font-family: monospace;'>QP</span><span>'), and mixed-integer quadratic programming ('</span><span style=' font-family: monospace;'>MIQP</span><span>').</span></div><h2  class = 'S6'><span style=' font-weight: bold;'>Model setup.</span></h2><div  class = 'S1'><span>Load the modelEcore, and define the uptake of nutrients by the modelEcore. The substrate used is glucose, and for this tutorial limit its uptake up to 18 mmol/(gDW·h).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'ecoli_core_model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >modelDirectory = getDistributedModelFolder(modelFileName); </span><span style="color: rgb(2, 128, 9);">%Look up the folder for the distributed Models.</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >modelFileName= [modelDirectory filesep modelFileName]; </span><span style="color: rgb(2, 128, 9);">% Get the full path. Necessary to be sure, that the right model is loaded</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >model = readCbModel(modelFileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >, -18, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><h2  class = 'S2'><span>PROCEDURE</span></h2><h2  class = 'S6'><span>Identify essentail reactions: perform a gene knocked-out analysis. </span></h2><div  class = 'S1'><span>If the modelEcore is not able to grow when a certain gene is knocked-out, the name of that gene will be saved as an essential gene. Even if a very small growth is calculated, the model will be considered as not growing and the gene will be recorded in an 'essential_genes' vector. Here no growth is defined as growth lower than 0.000001. The remaining non-essentail genes will be stored in a 'non_EG' vector.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >[grRatio, grRateKO, grRateWT, delRxns,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    hasEffect] = singleGeneDeletion(model, </span><span style="color: rgb(170, 4, 249);">'FBA'</span><span >, model.genes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >essential_genes = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >non_EG = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >tol = 1e-6;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >n = 1:length(grRateKO)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >(grRateKO(n)&lt;tol)||(isnan(grRateKO(n))==1)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        essential_genes = [essential_genes; model.genes(n)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        non_EG = [non_EG; n];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">% find essential reactions</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >RxnRatio = singleRxnDeletion(model);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >RxnRatio(isnan(RxnRatio)) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >pFBAEssentialRxns = model.rxns(RxnRatio &lt; tol);</span></span></div></div></div><h2  class = 'S2'><span>Identify non-essentail reactions that can or cannot carry flux:</span></h2><div  class = 'S1'><span>A FVA is performed without any biomass constraint. Therefore, for the </span><span style=' font-family: monospace;'>fluxVariability</span><span> function set the percentage of optimal solution to zero %. The reactions that do not carry flux will be stored in a vector called </span><span style=' font-family: monospace;'>pFBAnoFluxRxn</span><span> and the reaction that do carry flux in another vector called </span><span style=' font-family: monospace;'>pFBAfluxRxn</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >[minFluxglc, maxFluxglc] = fluxVariability(model, 0);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAnoFluxRxn = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAfluxRxn = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i=1:length(model.rxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >(abs(minFluxglc(i))&lt;tol)&amp;&amp;(abs(maxFluxglc(i))&lt;tol)  </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        pFBAnoFluxRxn = [pFBAnoFluxRxn i];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        pFBAfluxRxn = [pFBAfluxRxn i];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAfluxRxn = pFBAfluxRxn';</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >ZeroFluxRxns = model.rxns(pFBAnoFluxRxn) </span></span></div></div></div><div  class = 'S11'><span>Now, it is necessary to know which genes are associated to the reactions not carrying any flux. The </span><span style=' font-family: monospace;'>rxnGeneMat</span><span> filed in modelEcore stores information that connects genes to reactions. To extract information from </span><span style=' font-family: monospace;'>rxnGeneMat</span><span>, first converted it into a binary matrix of zeros and ones, afterwhich, use this matrix to get the names of the genes related with the </span><span style=' font-family: monospace;'>pFBAnoFluxRxn</span><span>. Then create a new vector, </span><span style=' font-family: monospace;'>pFBAfluxGenes</span><span>, of non essential genes that can carry flux and another vector, </span><span style=' font-family: monospace;'>pFBAnofluxGenes</span><span>, of non-essential genes that cannot carry flux.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >RxnGMat = full(model.rxnGeneMat);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAfluxGenes = non_EG;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAnoFluxGenes = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(pFBAnoFluxRxn)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    listGenes = find(RxnGMat(pFBAnoFluxRxn(i),:));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >n = 1:length(listGenes)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        pos = find(non_EG==listGenes(n));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >pos </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            pFBAnoFluxGenes = [pFBAnoFluxGenes; model.genes(non_EG(pos))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            pFBAfluxGenes(pos) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAnoFluxGenes = unique(pFBAnoFluxGenes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >pFBAfluxGenes(pFBAfluxGenes==0) = [];</span></span></div></div></div><h2  class = 'S2'><span>Identify MLE reactions:</span></h2><div  class = 'S1'><span>As suggested by Lewis et al.</span><span texencoding="$^1" style="vertical-align:-5px"><img src="" width="8" height="19" /></span><span>, calculate a FBA and set the FBA solution (i.e. optimal growth rate) as the lower bound of the objective function (in this case biomass production). The FBA is run with a </span><span>maxmial optimization of biomass production. Then set the optimal solution (</span><span style=' font-family: monospace;'>FBAsolution.f</span><span>) as the lower bound in the model using the </span><span style=' font-family: monospace;'>changeRxnBounds </span><span>function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >FBAsolution = optimizeCbModel(model, </span><span style="color: rgb(170, 4, 249);">'max'</span><span >); </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >model = changeRxnBounds(model, </span><span style="color: rgb(170, 4, 249);">'Biomass_Ecoli_core_w_GAM'</span><span >, FBAsolution.f, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><div  class = 'S11'><span>Then a FVA was run, but the percentage of optimal solution was set up to 95%. This simulation provides a minimum and a maximum flux balance solution that allows at least a 95% of the optimal solution for the objective function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: normal"><span >[minFlux2, maxFlux2] = fluxVariability(model,95);</span></span></div></div></div><div  class = 'S11'><span>The list of reactions carying flux will be scanned, and the ones that are "turned off" when the system is forced to achieve certain biomass production are MLE reactions. MLE reations will be stored in the vector, </span><span style=' font-family: monospace;'>RxnMLE</span><span>, and the remaining reactions will be stored in the vector, </span><span style=' font-family: monospace;'>restRxn</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >RxnMLE = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >restRxn = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(pFBAfluxRxn)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >(abs(minFlux2(pFBAfluxRxn(i)))&lt;tol)&amp;&amp;(abs(maxFlux2(pFBAfluxRxn(i)))&lt;tol);    </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        RxnMLE = [RxnMLE pFBAfluxRxn(i)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        restRxn = [restRxn pFBAfluxRxn(i)];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >RxnMLEname = model.rxns(RxnMLE)</span></span></div></div></div><h2  class = 'S2'><span>Identify Optimal and ELE reactions:</span></h2><div  class = 'S1'><span>Next run an FBA, calculating a minimal optimization </span><span>of biomass production. Then set the bounds of all reactions to it respective minimal flux balance solution (</span><span style=' font-family: monospace;'>FBAsolution.x</span><span>) using the </span><span style=' font-family: monospace;'>changeRxnBounds() </span><span>function. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >FBAsolution = optimizeCbModel(model,</span><span style="color: rgb(170, 4, 249);">'min'</span><span >,</span><span style="color: rgb(170, 4, 249);">'one'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >model = changeRxnBounds(model, model.rxns, FBAsolution.x, </span><span style="color: rgb(170, 4, 249);">'b'</span><span >);</span></span></div></div></div><div  class = 'S11'><span>Finally, run one last FVA for 100% of the optimal solution. The remaining</span><span> reactions in the </span><span style=' font-family: monospace;'>restRxn</span><span> variable were then clasified as Enzymatially Less Eficient Reactions (RxnELE), if the reactions cannot carry any flux, or as Optimal Reactions (RxnOptima), if they can carry flux.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >[minFlux3, maxFlux3] = fluxVariability(model, 100);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAopt_Rxns = model.rxns((abs(minFlux3)+abs(maxFlux3))&gt;=tol);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAopt_Rxns = unique(regexprep(pFBAopt_Rxns, </span><span style="color: rgb(170, 4, 249);">'_[f|b]$'</span><span >,</span><span style="color: rgb(170, 4, 249);">''</span><span >));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >pFBAopt_Rxns = setdiff(pFBAopt_Rxns, pFBAEssentialRxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >ELE_Rxns = model.rxns((abs(minFlux3)+abs(maxFlux3))&lt;=tol);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >ELE_Rxns = setdiff(ELE_Rxns, RxnMLEname);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >ELE_Rxns = setdiff(ELE_Rxns, ZeroFluxRxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >RxnELE = findRxnIDs(model, ELE_Rxns);</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >RxnOptima = findRxnIDs(model, pFBAopt_Rxns);</span></span></div></div></div><h2  class = 'S2'><span>Classify the genes:</span></h2><div  class = 'S1'><span>The last step is to associate the genes that are related with each reaction. The main point of this is to classify the genes into the 5 different groups (Figure 3) and store them into different vectors:</span></div><ol  class = 'S3'><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>Essential genes:</span><span> metabolic genes necessary for growth in the given media ('essential_genes').</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>pFBA optima: </span><span>non-essential genes contributing to the optimal growth rate and minimum gene-associated flux ('OptimaGenes').</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>Enzymatically less efficient (ELE): </span><span>genes requiring more flux through enzymatic steps than alternative pathways that meet the same predicted growth rate ('ELEGenes').</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>Metabolically less efficient (MLE):</span><span> genes requiring a growth rate reduction if used ('MLEGenes').</span></li><li  class = 'S4'><span style=' font-weight: bold; font-style: italic;'>pFBA no-flux:</span><span> genes that are unable to carry flux in the experimental conditions ('pFBAnoFluxGenes').</span></li></ol><div  class = 'S5'><img class = "imageNode" src = "" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S1'><span>                                                                                  Figure 3: Gene classes retrieved through pFBA.</span></div><div  class = 'S1'><span>Some genes may </span><span>not fit in any of this 5 categories. These genes will be saved in a vetor called 'remainingGenes'.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >OptimaGenes = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >restGenes = pFBAfluxGenes;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(RxnOptima)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    listGenes = find(RxnGMat(RxnOptima(i), :));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >n = 1:length(listGenes)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        pos = find(pFBAfluxGenes==listGenes(n));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >pos </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            OptimaGenes = [OptimaGenes; model.genes(pFBAfluxGenes(pos))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            restGenes(pos,1) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >OptimaGenes = unique(OptimaGenes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >restGenes(restGenes==0) = [];    </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >ELEGenes = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >restGenes2 = restGenes;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(RxnELE)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    listGenes = find(RxnGMat(RxnELE(i), :));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >n = 1:length(listGenes)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        pos = find(restGenes==listGenes(n));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >pos </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            ELEGenes = [ELEGenes; model.genes(restGenes(pos))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            restGenes2(pos, 1) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >ELEGenes = unique(ELEGenes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >restGenes2(restGenes2==0) = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >MLEGenes = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >finalRemainingGenes = restGenes2;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(RxnMLE)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    listGenes = find(RxnGMat(RxnMLE(i),:));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">for </span><span >n = 1:length(listGenes)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        pos = find(restGenes2==listGenes(n));</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">if </span><span >pos </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            MLEGenes = [MLEGenes; model.genes(restGenes2(pos))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >            finalRemainingGenes(pos, 1) = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        </span><span style="color: rgb(14, 0, 255);">end </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >MLEGenes = unique(MLEGenes);</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >finalRemainingGenes(finalRemainingGenes==0) = []; </span></span></div></div><div class="inlineWrapper"><div  class = 'S9'></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >remainingGenes = [];</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >n = 1:length(finalRemainingGenes)</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >        remainingGenes = [remainingGenes; model.genes(finalRemainingGenes(n))];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S6'><span style=' font-weight: bold;'>Pri</span><span style=' font-weight: bold;'>nt results:</span></h2><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S8'><span style="white-space: normal"><span >essential_genes</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >OptimaGenes</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >ELEGenes</span></span></div></div><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >MLEGenes</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >pFBAnoFluxGenes</span></span></div></div></div><h2  class = 'S2'><span>TIMING</span></h2><div  class = 'S1'><span>The tutorial runs in a few minutes.</span></div><h2  class = 'S2'><span>References</span></h2><div  class = 'S1'><span>[1] Lewis et al. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. </span><span style=' font-style: italic;'>Mol Syst Biol. </span><span>6:390 (2010).</span></div><div  class = 'S1'><span>[2] Orth, J., Fleming, R.M., Palsson B. Ø. Reconstruction and Use of Microbial Metabolic Networks: the Core </span><span>Escherichia coli</span><span> Metabolic Model as an Educational Guide. </span><span style=' font-style: italic;'>EcoSal Plus.</span><span> 4(1) (2010).</span></div><div  class = 'S1'><span>[3] Thiele, I., et al. A community-driven global reconstruction of human metabolism. </span><span style=' font-style: italic;'>Nat Biotechnol</span><span>. 31(5):419-425 (2013).</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% A step-by-step guide to parsimoneous enzyme usage Flux Balance Analysis - pFBA
% *Author(s): Francisco José Pardo Palacios, Ines Thiele, LCSB, University of 
% Luxembourg.*
% 
% *Reviewer(s): Sebastián Mendoza, Center for Mathematical Modeling, University 
% of Chile. Catherine Clancy, LCSB, University of Luxembourg.*
% 
% *Lin Wang (Costas D. Maranas Lab), Joshua Chan (Costas D. Maranas Lab), Chiam 
% Yu Ng (Costas D. Maranas Lab)*
%% INTRODUCTION
% This tutorial shows how the parsimoneous enzyme usage Flux Balance Analysis 
% (pFBA), as described in Lewis et al.$$^1$, has been implemented in The COBRA 
% Toolbox as the function |pFBA()|.
% 
% The main aim of the tutorial is to explain how the calculations are carried 
% out in order to understand the pFBA analysis, and to be able to classify, under 
% certain conditions, the genes of a model as: essential, pFBA optima, Enzymatically 
% Less Efficient (ELE), Metabolically Less Efficient (MLE) or pFBA no-flux genes 
% (Figure 1). 
%% 
% # _*Essential genes:*_ metabolic genes necessary for growth in the given media.
% # _*pFBA optima:*_ non-essential genes contributing to the optimal growth 
% rate and minimum gene-associated flux.
% # _*Enzymatically less efficient (ELE):*_ genes requiring more flux through 
% enzymatic steps than alternative pathways that meet the same predicted growth 
% rate.
% # _*Metabolically less efficient (MLE):*_ genes requiring a growth rate reduction 
% if used.
% # _*pFBA no-flux:*_ genes that are unable to carry flux in the experimental 
% conditions.
%% 
% 
% 
% Figure 1: Gene/enzyme classification scheme used by pFBA
% 
% This tutorial will use the _E. coli core_ reconstruction$$^2$ as the model 
% of choice, and will be called herein as |modelEcore|. The results obtained could 
% then be compared to data from evolved E. coli and observe if the |modelEcore| 
% can predict its evolution. In order to investigate this, all the steps described 
% in the pFBA flowchart (Figure 2) should be followed, and are demonstrated in 
% this tutorial. 
% 
% 
% 
% Figure 2: pFBA flowchart
%% EQUIPMENT SETUP
%% *Initialize the COBRA Toolbox.*
% If necessary, initialize The Cobra Toolbox using the |initCobraToolbox| function.

initCobraToolbox(false) % false, as we don't want to update
%% *Setting the* optimization *solver.*
% This tutorial will be run with a |'glpk'| package, which is a linear programming 
% ('|LP'|) solver. The |'glpk'| package does not require additional instalation 
% and configuration.

solverName = 'glpk';
solverType = 'LP'; 
changeCobraSolver(solverName, solverType);
% changeCobraSolver('ibm_cplex');
% changeCobraSolver('gurobi');
%% 
% However, for the analysis of larger models, such as Recon 2.04$$^3$, it is 
% not recommended to use the |'glpk'| package but rather an industrial strength 
% solver, such as the |'gurobi' or 'ibm_cplex'| package.
% 
% A solver package may offer different types of optimization programmes to solve 
% a problem. The above example used a LP optimization, other types of optimization 
% programmes include; mixed-integer linear programming ('|MILP|'), quadratic programming 
% ('|QP|'), and mixed-integer quadratic programming ('|MIQP|').
%% *Model setup.*
% Load the modelEcore, and define the uptake of nutrients by the modelEcore. 
% The substrate used is glucose, and for this tutorial limit its uptake up to 
% 18 mmol/(gDW·h).

modelFileName = 'ecoli_core_model.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
model = readCbModel(modelFileName);
model = changeRxnBounds(model, 'EX_glc(e)', -18, 'l');
%% PROCEDURE
%% Identify essentail reactions: perform a gene knocked-out analysis. 
% If the modelEcore is not able to grow when a certain gene is knocked-out, 
% the name of that gene will be saved as an essential gene. Even if a very small 
% growth is calculated, the model will be considered as not growing and the gene 
% will be recorded in an 'essential_genes' vector. Here no growth is defined as 
% growth lower than 0.000001. The remaining non-essentail genes will be stored 
% in a 'non_EG' vector.

[grRatio, grRateKO, grRateWT, delRxns,...
    hasEffect] = singleGeneDeletion(model, 'FBA', model.genes);
essential_genes = [];
non_EG = [];
tol = 1e-6;
for n = 1:length(grRateKO)
    if (grRateKO(n)<tol)||(isnan(grRateKO(n))==1)
        essential_genes = [essential_genes; model.genes(n)];
    else
        non_EG = [non_EG; n];
    end
end

% find essential reactions
RxnRatio = singleRxnDeletion(model);
RxnRatio(isnan(RxnRatio)) = 0;
pFBAEssentialRxns = model.rxns(RxnRatio < tol);
%% Identify non-essentail reactions that can or cannot carry flux:
% A FVA is performed without any biomass constraint. Therefore, for the |fluxVariability| 
% function set the percentage of optimal solution to zero %. The reactions that 
% do not carry flux will be stored in a vector called |pFBAnoFluxRxn| and the 
% reaction that do carry flux in another vector called |pFBAfluxRxn|.

[minFluxglc, maxFluxglc] = fluxVariability(model, 0);
pFBAnoFluxRxn = [];
pFBAfluxRxn = [];
for i=1:length(model.rxns)
    if (abs(minFluxglc(i))<tol)&&(abs(maxFluxglc(i))<tol)  
        pFBAnoFluxRxn = [pFBAnoFluxRxn i];
    else
        pFBAfluxRxn = [pFBAfluxRxn i];
    end
end
pFBAfluxRxn = pFBAfluxRxn';
ZeroFluxRxns = model.rxns(pFBAnoFluxRxn) 
%% 
% Now, it is necessary to know which genes are associated to the reactions not 
% carrying any flux. The |rxnGeneMat| filed in modelEcore stores information that 
% connects genes to reactions. To extract information from |rxnGeneMat|, first 
% converted it into a binary matrix of zeros and ones, afterwhich, use this matrix 
% to get the names of the genes related with the |pFBAnoFluxRxn|. Then create 
% a new vector, |pFBAfluxGenes|, of non essential genes that can carry flux and 
% another vector, |pFBAnofluxGenes|, of non-essential genes that cannot carry 
% flux.

RxnGMat = full(model.rxnGeneMat);
pFBAfluxGenes = non_EG;
pFBAnoFluxGenes = [];

for i = 1:length(pFBAnoFluxRxn)
    listGenes = find(RxnGMat(pFBAnoFluxRxn(i),:));
    for n = 1:length(listGenes)
        pos = find(non_EG==listGenes(n));
        if pos 
            pFBAnoFluxGenes = [pFBAnoFluxGenes; model.genes(non_EG(pos))];
            pFBAfluxGenes(pos) = 0;
        end 
    end
end

pFBAnoFluxGenes = unique(pFBAnoFluxGenes);
pFBAfluxGenes(pFBAfluxGenes==0) = [];
%% Identify MLE reactions:
% As suggested by Lewis et al.$$^1$, calculate a FBA and set the FBA solution 
% (i.e. optimal growth rate) as the lower bound of the objective function (in 
% this case biomass production). The FBA is run with a maxmial optimization of 
% biomass production. Then set the optimal solution (|FBAsolution.f|) as the lower 
% bound in the model using the |changeRxnBounds| function.

FBAsolution = optimizeCbModel(model, 'max'); 
model = changeRxnBounds(model, 'Biomass_Ecoli_core_w_GAM', FBAsolution.f, 'l');
%% 
% Then a FVA was run, but the percentage of optimal solution was set up to 95%. 
% This simulation provides a minimum and a maximum flux balance solution that 
% allows at least a 95% of the optimal solution for the objective function.

[minFlux2, maxFlux2] = fluxVariability(model,95);
%% 
% The list of reactions carying flux will be scanned, and the ones that are 
% "turned off" when the system is forced to achieve certain biomass production 
% are MLE reactions. MLE reations will be stored in the vector, |RxnMLE|, and 
% the remaining reactions will be stored in the vector, |restRxn|.

RxnMLE = [];
restRxn = [];
for i = 1:length(pFBAfluxRxn)
    if (abs(minFlux2(pFBAfluxRxn(i)))<tol)&&(abs(maxFlux2(pFBAfluxRxn(i)))<tol);    
        RxnMLE = [RxnMLE pFBAfluxRxn(i)];
    else
        restRxn = [restRxn pFBAfluxRxn(i)];
    end
end

RxnMLEname = model.rxns(RxnMLE)
%% Identify Optimal and ELE reactions:
% Next run an FBA, calculating a minimal optimization of biomass production. 
% Then set the bounds of all reactions to it respective minimal flux balance solution 
% (|FBAsolution.x|) using the |changeRxnBounds()| function. 

FBAsolution = optimizeCbModel(model,'min','one');
model = changeRxnBounds(model, model.rxns, FBAsolution.x, 'b');
%% 
% Finally, run one last FVA for 100% of the optimal solution. The remaining 
% reactions in the |restRxn| variable were then clasified as Enzymatially Less 
% Eficient Reactions (RxnELE), if the reactions cannot carry any flux, or as Optimal 
% Reactions (RxnOptima), if they can carry flux.

[minFlux3, maxFlux3] = fluxVariability(model, 100);

pFBAopt_Rxns = model.rxns((abs(minFlux3)+abs(maxFlux3))>=tol);
pFBAopt_Rxns = unique(regexprep(pFBAopt_Rxns, '_[f|b]$',''));
pFBAopt_Rxns = setdiff(pFBAopt_Rxns, pFBAEssentialRxns)
ELE_Rxns = model.rxns((abs(minFlux3)+abs(maxFlux3))<=tol);
ELE_Rxns = setdiff(ELE_Rxns, RxnMLEname);
ELE_Rxns = setdiff(ELE_Rxns, ZeroFluxRxns)
RxnELE = findRxnIDs(model, ELE_Rxns);
RxnOptima = findRxnIDs(model, pFBAopt_Rxns);
%% Classify the genes:
% The last step is to associate the genes that are related with each reaction. 
% The main point of this is to classify the genes into the 5 different groups 
% (Figure 3) and store them into different vectors:
%% 
% # _*Essential genes:*_ metabolic genes necessary for growth in the given media 
% ('essential_genes').
% # _*pFBA optima:*_ non-essential genes contributing to the optimal growth 
% rate and minimum gene-associated flux ('OptimaGenes').
% # _*Enzymatically less efficient (ELE):*_ genes requiring more flux through 
% enzymatic steps than alternative pathways that meet the same predicted growth 
% rate ('ELEGenes').
% # _*Metabolically less efficient (MLE):*_ genes requiring a growth rate reduction 
% if used ('MLEGenes').
% # _*pFBA no-flux:*_ genes that are unable to carry flux in the experimental 
% conditions ('pFBAnoFluxGenes').
%% 
% 
% 
% Figure 3: Gene classes retrieved through pFBA.
% 
% Some genes may not fit in any of this 5 categories. These genes will be saved 
% in a vetor called 'remainingGenes'.

OptimaGenes = [];
restGenes = pFBAfluxGenes;
for i = 1:length(RxnOptima)
    listGenes = find(RxnGMat(RxnOptima(i), :));
    for n = 1:length(listGenes)
        pos = find(pFBAfluxGenes==listGenes(n));
        if pos 
            OptimaGenes = [OptimaGenes; model.genes(pFBAfluxGenes(pos))];
            restGenes(pos,1) = 0;
        end 
    end
end
OptimaGenes = unique(OptimaGenes);
restGenes(restGenes==0) = [];    

ELEGenes = [];
restGenes2 = restGenes;
for i = 1:length(RxnELE)
    listGenes = find(RxnGMat(RxnELE(i), :));
    for n = 1:length(listGenes)
        pos = find(restGenes==listGenes(n));
        if pos 
            ELEGenes = [ELEGenes; model.genes(restGenes(pos))];
            restGenes2(pos, 1) = 0;
        end 
    end
end
ELEGenes = unique(ELEGenes);
restGenes2(restGenes2==0) = [];

MLEGenes = [];
finalRemainingGenes = restGenes2;
for i = 1:length(RxnMLE)
    listGenes = find(RxnGMat(RxnMLE(i),:));
    for n = 1:length(listGenes)
        pos = find(restGenes2==listGenes(n));
        if pos 
            MLEGenes = [MLEGenes; model.genes(restGenes2(pos))];
            finalRemainingGenes(pos, 1) = 0;
        end 
    end
end
MLEGenes = unique(MLEGenes);
finalRemainingGenes(finalRemainingGenes==0) = []; 

remainingGenes = [];
for n = 1:length(finalRemainingGenes)
        remainingGenes = [remainingGenes; model.genes(finalRemainingGenes(n))];
end
%% *Print results:*

essential_genes
OptimaGenes
ELEGenes
MLEGenes
pFBAnoFluxGenes
%% TIMING
% The tutorial runs in a few minutes.
%% References
% [1] Lewis et al. Omic data from evolved E. coli are consistent with computed 
% optimal growth from genome-scale models. _Mol Syst Biol._ 6:390 (2010).
% 
% [2] Orth, J., Fleming, R.M., Palsson B. Ø. Reconstruction and Use of Microbial 
% Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational 
% Guide. _EcoSal Plus._ 4(1) (2010).
% 
% [3] Thiele, I., et al. A community-driven global reconstruction of human metabolism. 
% _Nat Biotechnol_. 31(5):419-425 (2013).
##### SOURCE END #####
-->
</div></body></html>